Method for detecting and discriminating breathing patterns from respiratory signals

ABSTRACT

A signal representative of a patient&#39;s respiration is split into equal length epochs. A primary feature is extracted from each epoch that acts as a compressed representation of the signal events. Statistics are applied to the primary feature to produce one or more secondary features that represent the entire epoch. Each secondary feature is grouped with one or more other features that are extracted from the entire epoch rather than selected epoch events. This grouping is the epoch pattern. The pattern is manipulated with suitable classifier algorithm to produce a probability for each class within the algorithm, that the signal may be representative of a disease state (Cheyne-Stokes, OSA etc). The epoch is assigned to the class with the highest probability. Also defined are methods of detecting Cheyne-Stokes breathing by analyzing a signal to detect one or regions of hyperpnoea and if the length of a hyperpnoea exceeds a parameter, Cheyne-Stokes is present.

CROSS-REFERENCES TO RELATED APPLICATIONS

This application claims the priority of U.S. Provisional PatentApplication 60/638,169, filed Dec. 23, 2004.

FIELD OF THE INVENTION

The present invention relates to methods, algorithms and diagnosticapparatus for

the detection of sleep disordered breathing patterns and thediscrimination between patterns of different disease states such asobstructive sleep apnea, central sleep apnea and Cheyne-Stokes breathingand mixed sleep apnea.

BACKGROUND OF THE INVENTION

Breathing Disorders

Sleep-disordered breathing (SDB) encompasses a group of disorders wherethe breathing pattern or quality of ventilation is abnormal duringsleep. Obstructive sleep apnea (OSA), the most common such disorder(effecting possible 4-5% of the adult population), is characterized byrepetitive closing or collapse of the upper airway and partial orcomplete diminution of breathing. The obstruction is normally ended bythe patient arousing briefly when the muscles of the upper airway act toclear the obstruction. During the repetitive cycle of obstruction andarousal, the OSA patient will always continue to make “efforts” tobreath; in other words there is no central or brain-mediated disruptionto the breathing cycle.

Conversely, in central sleep apnea (CSA), there is a disruption tobreathing which is brain or control-centre in origin. Cheyne-Stokesbreathing (CS) is one of the more common forms of CSA. It is caused byan abnormal limit cycle instability of the patient's respiratorycontroller in which there are rhythmic alternating periods of waxing andwaning ventilation. Patients with cardiac failure (a condition where theheart fails to pump adequately) often have CSA, especially as thecondition deteriorates or where therapy has ceased to allow effectivecompensation by the heart. Cheyne-Stokes breathing appears as a cyclicalvariation in tidal volume seen in heart failure patients. The cycleconsists of an apnoea or hypopnoea followed by an overshootinghyperpnoea which often (but not always) has a characteristic hump backedmorphology s a.k.a. a “Sydney Harbor Bridge” shape. The exact cause ofCS breathing is not fully understood. However, the characteristic waxingand waning cycle is strongly reminiscent of limit cycles in a poorlyadjusted control system with a maladjusted gain or destabilizingfeedback-loop delay.

Sleep-disordered breathing is undesirable in all its forms because itdisrupts sleep architecture (the pattern and proportion of the differentforms of sleep) leading to daytime somnolence. The repetitive cessationor diminution of ventilation causes (sometimes dramatic) drops in bloodoxygenation levels. These and other complications are probablyresponsible for the now established sequelae of cardiovascularconditions.

The treatment of choice for OSA is continuous positive airway pressure(CPAP) as first described by Sullivan [Sullivan C E, et al. Reversal ofobstructive sleep apnea by continuous positive airway pressure appliedthrough the nares. Lancet Apr. 18, 1981 ;1(8225):862-5]. CPAP is alsoused to treat some heart-failure patients with CSA and congestive heartfailure (fluid on the lungs). However, Cheyne-Stokes breathing isineffectively treated by CPAP and may require the application ofservo-ventilation [Teschler H et al. Adaptive pressure supportservo-ventilation: a novel treatment for Cheyne-Stokes respiration inheart failure. Am J Respir Crit Care Med. Aug. 15, 2001; 164(4):614-9.Berthon-Jones Ventilatory assistance for treatment of cardiac failureand Cheyne-Stokes breathing. U.S. Pat. No. 6,532,959].

Diagnosis from Multiple Signals

The gold standard for the diagnosis of SDB and sleep apnea is thepolysomnograph (PSG): the measurement and recording of a multitude ofphysiological signals during a stay overnight in a sleep laboratory.Briefly, the PSG signal ensemble normally includes one or more signalsindicative of a respiratory parameters such as patient airflow rate (forthe calculation of ventilation and the detection of apneas andhypopnoeas), multiple electroencephalogram (EEG), electrooculogram (EOG)and electromyogram (EMG) signals (for the determination of patient sleepstate, position and the detection of arousals from sleep), breathingeffort signals (either chest and abdominal distension bands or anesophageal pressure-measuring catheter), snore amplitude, and oxygensaturation. Another method of diagnosing SDB is polygraphy (PG) wherebya reduced number of parameters are recorded while the patient sleeps.These parameters include: nasal/oral airflow rate, snore amplitude,oxygen saturation, respiratory effort (thoracic and abdominal) and bodyposition.

In both the PSG and the PG a breathing-effort signal is recorded toenable the discrimination of OSA events from CSA or Cheyne-Stokesbreathing. (A third type of event is also possible-the mixed apnea-wherethe event is initiated by a centrally-mediated lack of breathing driveand ends with an airway obstruction and subsequent arousal). It isimpossible for the inexperienced observer to reliably determine the typeof apnea without reference to at least the flow signal and a measure ofbreathing effort. However, an experienced and trained observer (expert)can often readily detect patterns in a run of events (apneas/hypopnoeas)allowing a reliable determination of the type of underlying disease.This is especially true of Cheyne-Stokes breathing which has a verycharacteristic waxing and waning pattern of ventilation.

Simple Recording Devices

The performance of either a PSG or PG requires trained technicians, isexpensive, is time consuming and can itself introduce sleepdisturbances. Also, it is well known that a shortage of sleeplaboratories is hampering the diagnosis and treatment of current SDBpatients, let alone what is considered a vast undiagnosed population.For these reasons a type of “screening” device (e.g., the MICROMESAN®from MAP of Germany, or the APNEALINK™ from ResMed) is available to testpatients suspected of having sleep-disordered breathing. Such devicesare small, recording just one or two physiological signals, and can bereadily sent home with the patient for a screening study. For example:patients' nasal airflow can be recorded and later examined by aphysician using a personal computer and a connection to the device. Asoftware package would then be available to read the data from thedevice, show statistics and make recommendations regarding suspectedsleep-related pathology.

Diagnosis Classifier

The calculation of the apnea-hypopnoea index (AHI, number of such eventsper hour on average) is a measure regularly used to guide the directionof either treatment or further investigation with a fill PSG or PG. Acomputer program or algorithm which further enables the discriminationbetween different underlying disease states based on the recordedbreathing patterns provides added guidance to the clinical pathway. Astrong indication of Cheyne-Stokes disease, for example, would suggestcompletely different follow-up compared to the more common forms ofsleep apnea.

The concept of a classifier is common to many fields where it isdesirable to assign an object or an underlying state of an object to oneof a number of classes. This concept is used, for example, in the fieldsof voice recognition (where sound bytes are classified as differentwords or syllables), radar detection (where visual signals areclassified as enemy/friendly targets) and medical diagnosis (where testresults are used to classify a patient's disease state). The design of aclassifier falls under the field of Pattern Recognition and a classifiercan be of the supervised type (the classifier is built from trainingdata which has been pre-classed by a supervisor or “expert”) orunsupervised type (where the natural ordering or clustering of the datadetermines the different classes). Time signal classification usuallyrelies on representing the signal at particular time points with“features”. Features are simply numbers that distill the essence of thesignal at a point in time, a form of compression. A set (or vector) offeatures is called a “pattern”. A classifier takes a pattern andmanipulates it mathematically with a suitable algorithm to produce aprobability value for each of a number of classes. The pattern isassigned to the class with the highest probability.

In U.S. Pat. No. 6,839,581 there is disclosed a method for detecting CSrespiration in patients with congestive heart failure by performingspectral analysis of overnight oximeter recordings to obtain a set ofparameters that can be used in the construction of a classification treeand a trained neural network.

In summary, sleep-disordered breathing is a common syndrome withdifferent underlying disease types requiring very different treatmentoptions. There is a need for a small and relatively inexpensivescreening devices that can help unblock the treatment bottleneck thatcurrently exists at the sleep laboratory. An algorithm and diagnosticapparatus that can replicate the expert's ability to detect breathingpatterns associated with particular disease states will enhance thediagnosis and treatment of patients being screened for sleep-disorderedbreathing, or for monitoring patients already undergoing therapy. Whatis needed is an algorithm for flow data in the form of classifier.

What is particularly desirable is a method and apparatus for diagnosingCheyne-Stokes breathing from flow readings or oximeter readings by useof appropriate software in conjunction with a small hand-held device foruse in a home setting.

BRIEF SUMMARY OF THE INVENTION

The CS diagnosis system of the present invention uses patternclassification techniques on a digital computer to identify periods ofCS-like breathing by examining the flow signal alone. Ordinarily thedefinitive diagnosis of CS breathing relies on an “effort” signal,either esophageal pressure or an elastic band signal from the abdomen orthorax. An absence of effort denotes a central apnoea which mayotherwise be difficult to distinguish from an obstructive apnoea or amixed apnoea. A mixed apnoea is comprised of a central beginning(without effort) followed by a section of obstructed breaths once drivereturns.

APNEALINK™ nasal flow data without other channels is processed toclassify it as unambiguously Cheyne-Stokes (CS) breathing or nearly soand to then display a likely record to the physician for quick expertconfirmation. An APNEALINK™ recorder is a single channel battery-poweredrespiratory pressure sensor system and provides recordings ofrespiratory pressure during sleep. The APNEALINK™ is a small (hand held)device manufactured by ResMed, designed for use in a home setting whereit is worn strapped to the patient's chest. The device only recordsnasal flow (indirectly) using a nasal pressure-sensing catheter. Allrelevant respiratory information during sleep will be collected vianasal pressure cannula. This will allow cardiologists to manage suchpatients more expediently. For example, CS patients would go on to afull polysomnogram (PSG) workup for possible AutoSet CS therapy asappropriate. Non-CS patients might just go on to AutoSet to treat theunderlying OSA as appropriate. After suitable offline processing of thenasal flow signal using a PC, the following events can be detected anddisplayed: apnoeas, hypopnoeas, flow-limitation and snore.

Cheyne-Stokes Detection Algorithm

The CS-detection algorithm uses the nasal flow signal from a device suchas ResMed's APNEALINK® or other signal indicative of at least onerespiratory parameter together with pattern recognition techniques toassign a probability of CS breathing to each 30 minute epoch of flowrecorded. This invention details the initial filtering and “event”detection, where events are defined as regions of hypopnoea-hyperpnoeasequence characteristic of CS breathing. The detection of such eventsmay be determined from the duration of one or more regions of hyperpnoeawhen the duration of the hyperpnoea exceeds a threshold or a statisticof the duration of regions of hyperpnoea exceeds a threshold. Such astatistic may be an average or a standard deviation or other statisticsspecified below.

Pattern classification techniques are statistical and rely on aso-called “training” data set by which a “classifier” can be trained torecognize certain “patterns”, in this case CS breathing. A pattern is agroup or vector of features. A feature is a number which represents someaspect of the signal being examined. An example of a feature is apnoealength. A pattern might be the group comprising apnoea length,hyperpnoea length and a number representing the closeness of the shapeof the hyperpnoea to a “harbor bridge”.

One aspect of the invention is directed to a method and apparatus orsystem capable of better diagnosing the presence of sleep disorders,preferably with a higher level of confidence. The diagnosis may compriseanalyzing a signal indicative of a respiratory parameter to determine arate of increase of the signal in the region from hypopnea to hyperpnoeaand where the rate of increase is a slow increase, concluding thatCheyne-Stokes breathing is present and where the rate of increase is asudden increase, concluding that Cheyne-Stokes breathing is absent.

According to one aspect of the invention, a signal representative of apatient's respiration is split into equal length epochs which can be aslong (the entire record) or as short (the length of a representativehypopnoea-hyperpnoea sequence) as desired. Preferably, the signal willbe subject to a number of pre-processing steps in order to filter outnoise and zero the baseline, for example.

Preferably, from each epoch one or more primary features is extractedfrom the signal that act as a compressed representation of the signalevents. By events it is meant: e.g., apneas, hypopnoeas and hyperpnoeas.Statistics are applied to the primary feature(s) to produce one or moresecondary features which represent the entire epoch. Each secondaryfeature is grouped with one or more other features that is extractedfrom the entire epoch rather than from the epoch events. This finalgroup of features is the epoch pattern.

The epoch pattern is preferably manipulated with a suitable classifieralgorithm to produce a probability for each possible class that thesignal may be representative of (e.g. Cheyne-Stokes breathing, OSAetc.). The epoch is assigned to the class with the highest probabilityand the class and the strength of the probability can be reported as anindication of the underlying disease state.

The classifier algorithm is preferably learned from a training data setwhich has been pre-classified by a human expert. In this sense theclassifier is of the supervised machine learning type.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of the signal processing pathway leading fromthe patient's respiratory signal through pre-processing, featureextraction based on epochs, through to classification.

FIG. 2 shows a typical respiratory signal epoch including a number of“events” (in this case apnea-hyperpnoea sequences). Several primaryfeatures are either shown explicitly (2.1 apnea/hypopnoea length, 2.2hyperpnoea length, 2.3 cycle length) or inferred (2.4 shape features ofthe hyperpnoea, 2.5 a feature representing the initial “jump” at thebeginning of the hyperpnoea).

FIG. 3 shows details of the calculation of primary features. 3.1 showsthe basis fictions used in the determination of the hyperpnoea shapefeatures. 3.2 shows an hyperpnoea typical of OSA together with thecalculation of the jump feature displayed graphically. 3.3 shows asimilar depiction of an hyperpnoea more typical of CS breathing. In bothcases the calculated jump features and shape features are tabled.

FIG. 4 shows examples of epoch classification, e.g., using bar charts.

FIG. 5 shows the distribution of the normalized max jump in a hyperpnoeasignal.

FIG. 6 shows a cluster analysis of CS and OSA patients.

FIG. 7 shows results from a LD (linear) classifier.

FIG. 8 shows results from a QD (quadratic) classifier.

FIG. 9 shows the correction of data for baseline offset.

FIG. 10 shows a Cheyne-Stokes flow waveform, the long-term ventilationand the left-shifted 10-second ventilation for a typical patient.

FIG. 11 depicts those parts of the flow waveform that correspond tohyperpnoeas.

FIG. 12 shows a Cheyne-Stokes patient's nasal flow signal over about 15minutes.

FIG. 13 shows a patient's SpO2 signal (saturation) and ventilationsignal (low-pass filtered absolute value of flow).

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

Process Description

FIG. 1 shows one embodiment of the classification process. While thefollowing may be explained in terms of a sequential process; it isunderstood that the process can be carried out using a non-linear,non-sequential, or non-staged process, or the order of the process maybe changed. Also, while FIG. 1 describes an entire process, aspects ofthe invention may relate to only a subset of that process. A signalrepresentative of respiration is first recorded from a patient using alogging device which includes a data-acquisition system and a memory.The respiratory signal is then processed either on-board by therecording device or off-line using a computer.

Preferably, the signal is initially pre-processed. For example, thesignal is filtered to remove unwanted noise and, where appropriate, thebaseline is zeroed. The signal may also be linearised depending on thetransducer used to detect the respiration.

In the next stage the signal is divided into n epochs of equal length.The epoch length can be as long as the entire record or as short as ispracticable to enable detection of respiratory patterns. In onepreferred embodiment the epoch length is 30 minutes.

FIG. 2 shows a typical epoch recording in a patient with Cheyne-Stokesbreathing. The shape of the curve is reminiscent of the shape of theSydney Harbor Bridge and is sometimes referred to by that name. See alsoFIG. 10. The recording consists of five “events”, each event consistingof a hypopnoea (in this case also an apnea) followed by a hyperpnoea.For each event an algorithm is used to detect the beginning and endpoints such that event lengths can be calculated: e.g., apnea/hypopnoealength and hyperpnoea length. A further algorithm may be used to rejectevents if they do not follow the correct sequence of hypopnoea/apneahyperpnoea. Another further algorithm may be used to reject events thatfall outside sensible length scale limits.

Determination of Shape Features

Each hyperpnoea is further processed to derive four so-called “shapefeatures”. These features indicate different shaped hyperpnoeas(bell-shaped versus triangle-shaped for example). The shape features arecalculated using singular value decomposition of the hyperpnoeaventilation signal as follows: First, the hyperpnoea is extracted fromthe respiratory signal and the absolute value is taken of therespiratory signal, giving a ventilation signal. The ventilation signalis scaled by its mean value to give a vector of values V_(hyperp). Formathematical convenience the time base of the hyperpnoea [0 . . . T],where T is the end of the hyperpnoea, is mapped to the interval [0 . . .2π]. A set of four orthogonal functions are calculated and arranged as a4×m matrix (where m is the number of values in the hyperpnoea signal). Aconvenient set of orthonormal function are:

$M_{Basis} = \begin{pmatrix}{\frac{1}{\sqrt{\pi}}{\sin\left( \frac{t}{2} \right)}} \\{\frac{1}{\sqrt{\pi}}{\cos\left( \frac{t}{2} \right)}} \\\frac{\left( {{3\pi\;{\sin(t)}} - {8\;{\cos\left( \frac{t}{2} \right)}}} \right)}{\sqrt{\pi\;\left( {{9\pi^{2}} - 64} \right)}} \\\frac{\left( {{3\pi\;{\cos(t)}} + {4\;{\sin\left( \frac{t}{2} \right)}}} \right)}{\sqrt{\pi\;\left( {{9\pi^{2}} - 16} \right)}}\end{pmatrix}$where t is the time base over the hyperpnoea from 0 to 2π. The basisfunctions are shown plotted in figure 3.1. The four shape features arethen calculated as:F _(P(1-4)) =V _(hyperp)×PseudoInverse(M _(Basis)),and are normalized by:

$F_{P{({1 - 4})}} = \frac{F_{P{({1 - 4})}}}{L_{2}{F_{P{({1 - 4})}}}}$where L₂∥ is the L2 or Euclidean norm.,

$\sqrt{\sum\limits_{i}^{n}\; x_{i}^{2}}$

The pseudoinverse M⁺ of a matrix M is a generalization of a matrixinverse, and exists for any (m,n) matrix M (for convenience assume m>n).If such a matrix M has full rank (n) one defines: M⁺=(M^(T)M)⁻¹M^(T).The solution of Mx=b is then x=M⁺b. (Pseudoinverses are useful becauseof a general theorem stating that F=M⁺v is the shortest length leastsquares solution to the problem MF=v.)

Jump Determination

Since sudden jumps in the ventilation/flow at the beginning of anhyperpnoea are characteristic of OSA, (see FIG. 2) each hyperpnoea isfurther processed to derive the so-called “jump” feature, indicative ofthe extent of any sudden increase in flow, as follows: Again, thehyperpnoea is extracted from the respiratory signal, the absolute valueis taken of the respiratory signal, giving a ventilation signal, adroopy peak-detector is used to approximate the envelope of theventilation signal:

e[1] = v[1] for  i = 2…  m if  v[i] ≥ e[i − 1] e[i] = v[i] else${e\lbrack i\rbrack} = {{e\left\lbrack {i - 1} \right\rbrack} + {\frac{1}{2.5f_{s}}\left( {{v\lbrack i\rbrack} - {e\left\lbrack {i - 1} \right\rbrack}} \right)}}$endwhere e[i] is the approximate envelope, f_(s) is the sampling frequencyand <[i] is the ventilation signal. The envelope is interpolated over anew two-second time base (chosen to be roughly the time-length of abreath) to give e₁ (between non-breathing intervals). The maximumpositive difference e_(1i)−e_(1(i−1)) (over the two second interval) isfound in the interpolated signal in the interval between the beginningof the envelope and the point at which the envelope attains its maximumvalue. Finally, the maximum difference is scaled by the mean value ofthe ventilation signal to give the “jump feature”. Figures 3.2 and 3.3show this process graphically for two representative hyperpnoeas.Secondary Feature Determination

Secondary features are calculated from primary features using the(measure of variation) statistics detailed below. (Note log denotes thelogarithm to base e.) First we define the standard deviation as:

${{{STD}\left( F_{P} \right)} = \sqrt{\frac{1}{n - 1}{\sum\limits_{i = 1}^{n}\;\left( {F_{Pi} - \overset{\_}{F_{P}}} \right)^{2}}}},{{{where}\mspace{14mu}\overset{\_}{F_{P}}} = {\frac{1}{n}{\sum\limits_{i = 1}^{n}\; F_{Pi}}}}$For length measures (e.g. hypopnoea length) and the jump feature thefour features are:

$1.\mspace{20mu}\frac{1}{n}{\sum\limits_{i = 1}^{n}\;{\log\left( F_{Pi} \right)}}$2.  log (STD(log (F_(P))))$3.\mspace{20mu}{\log\left( \sqrt{\frac{1}{n}{\sum\limits_{i = 1}^{n}\;\left( {F_{Pi} - F_{P{({i - 1})}}} \right)^{2}}} \right)}\left( {{To}\mspace{14mu}{get}\mspace{14mu} a\mspace{14mu}{normed}\mspace{14mu}{deviation}} \right)$4.  log (STD(F_(Pi) − F_(P(i − 1))))

For hyperpnoea shape features the four features are:

$1.\mspace{20mu}\frac{1}{n}{\sum\limits_{i = 1}^{n}\;\left( F_{Pi} \right)}$2.  log (STD((F_(P))))$3.\mspace{14mu}{\log\left( \sqrt{\frac{1}{n}{\sum\limits_{i = 1}^{n}\;\left( {F_{Pi} - F_{P{({i - 1})}}} \right)^{2}}} \right)}$4.  log (STD(F_(Pi) − F_(P(i − 1))))Additional Feature Determination

Additional features can be calculated using the entire (e.g. 30 minute)epoch signal. One such feature is derived from the spectrogram of theepoch signal and determining that Cheyne-Stokes breathing is present ifthe spectrogram indicates that the signal has a peak. This feature iscalculated as follows: First, the mean of the respiratory signal iscalculated and subtracted from the respiratory signal and the resultingsignal is chopped into n slices which overlap each other by exactly halfthe slice length. Each slice is next windowed, preferably using aHanning window (to reduce edge effects).

The use of a Hanning window to prepare the data for a FFT is as follows:The FFT function treats the N samples that it receives as though theyformed the basic unit of a repetitive waveform: It assumes that if onetook more samples they would repeat exactly, with the (N+1) sample beingidentical to the first sample, and so on. The usual case is that ifone's N samples start at one point in a cycle, they end at some otherpoint, so that if one really did play these back in a loop one would geta discontinuity between the last sample and the first. Hence one tapersboth ends of the set of samples down to zero, so they always line upperfectly if looped. The formal name for this process is “windowing”,and the “window function” is the shape that we multiply the data by.When the window function is the “raised cosine” 1+cos t the window istermed a Hanning window. Other periodic functions can be used, yieldingother windows.

Next, since CS data appears periodic, a fast Fourier transform isapplied to each windowed slice, yielding a complex vector result foreach slice. The absolute value is taken of each complex result yieldinga real valued vector per slice. The mean is taken of the resultingvectors to yield one vector. The natural log is taken of the subsequentvector and the values in the frequency range 0 Hz to 0.075 Hz areextracted to form a sub-vector, which is then de-trended. Cheyne-Stokesbehavior is present if the spectrogram indicates the signal has a peakin the range 0 Hz to 0.075 Hz.

Briefly, the method of detrended fluctuation analysis is useful inrevealing the extent of long-range correlations in time series, wherethe time series is a vector of data pairs (t_(i), x_(i)), where trepresents time and x represents the variable being measured.De-trending consists of subtracting from the x values, values that havebeen calculated using a polynomial of order n that has been fitted tothe data. For example, for order zero the polynomial is simply the meanof all the x values, and that mean is subtracted from the originalvalues. For order one, the polynomial is simply a linear fit to the data(t_(i), x_(i)). Values calculated using the best linear fit are thensubtracted from the original values (so removing any linear “trend”).For order two the fitted polynomial is a quadratic, for order three acubic etc.

The feature is then calculated as the maximum minus the mean of thede-trended vector. Alternatively one could calculate the entropy of theFFT instead of its peak.

Additional features can be derived by applying wavelet analysis to eachepoch. In this case wavelet coefficients or statistics derived fromwavelet coefficients are used as features for the epoch. This yields thelocation of the peak in time. In wavelet analysis a wave packet, offinite duration and with a specific frequency, is used as a windowfunction for an analysis of variance. This “wavelet” has the advantageof incorporating a wave of a certain period, as well as being finite inextent. A suitable wavelet (called the Morlet wavelet) is a sine wavemultiplied by a Gaussian envelope.

Classification

A subset of features is then selected for use by the classifier. It isknown that a particular subset of features can provide more accurateclassification than the full set of features. This is caused in part bythe so-called “curse of dimensionality”, whereby the required number oftraining samples increases with the number of features used. The curseof dimensionality causes networks with lots of irrelevant inputs tobehave relatively badly: Where the dimension of the input space is high,the network uses almost all its resources to represent irrelevantportions of the space.

An algorithm is employed to select the best subset based on the trainingdata. Ideally every subset of features should be tested for accuracy andthe best subset chosen. The number of subsets is 2^(n)−1 where n is thenumber of features. Unless there is a small number of features theexploration of all subsets is impractical and, in any case, accuracymeasures tend to be noisy which further hampers the search for the bestsubset. Alternative algorithms that enable selection of “good” featuresubsets include “best first”, “remove worst”, “random start with add andremove”, “simulated annealing” and genetic algorithms.

A method often used to measure accuracy is 10-fold cross-validation. Thetraining data are split into ten groups or folds and ten differentaccuracy tests are performed. In each case 9 tenths of the folds areused for training and the resulting classifier is tested for accuracy onthe remaining tenth. Statistics are performed on the 10 results to givea measure of accuracy.

Training the Classifier

Once a feature subset is chosen, the classifier is trained using theentire training data set. A number of classifier types are availableincluding: Baysean maximum likelihood linear and quadraticdiscriminants, neural networks and support vector machines. In each casea discriminant function is calculated which, when applied to featurescalculated from new data to be classified, provides probabilityestimates for different classes. The data (epoch) is assigned to theclass with the highest probability.

In one particular embodiment the discriminant function includes orpreferably consists of two weight vectors (of the same length as thefeature subset) and two constants. When the desired feature subset hasbeen extracted from the respiratory epoch, the discriminant functionsand probability are calculated as follows:

d₁ = W₁• F + C₁ d₂ = W₂• F + C₂${probability} = \frac{{\mathbb{e}}^{({{d\; 1} - {d\; 2}})}}{1 + {\mathbb{e}}^{({{d\; 1} - {d\; 2}})}}$

-   -   where W₁, W₂ are vectors and C₁, C₂ are constants.

The probability cutoff may be set at 0.5 in which case a probability of1.0 would equate to class A and a probability of 0.0 to class B. Thecutoff can be adjusted to suit the desired sensitivity and specificity.This is a two-way classification. With suitable training data, athree-way classification is also possible as are even higher n-wayclassifications.

In one particular embodiment the classification of each epoch could bedisplayed in a bar chart as in FIG. 4. Frame 4.1 shows a bar chart wheremany epochs show a high probability of a class of respiration (in thiscase CS-like breathing). This provides an “at-a-glance” indication of apatient record. Frame 4.2 shows a bar chart where only a single epochdisplays strong CS-like tendency. This provides an indication of wherein the patient's record a more detailed investigation is warranted.

Cheyne-Stokes Classifier Based on a Flow Signal or an Spo2 Signal orBoth

The APNEALINK device is capable of measuring an estimate of a patient'sflow signal which can be used as an input to the algorithm describedherein. Equally there are similar portable devices that can measure andlog SpO2, the saturation of oxyhemoglobin in the blood as estimated bypulse oximetry. Pulse oximetry is a simple non-invasive method ofmonitoring the percentage of haemoglobin (Hb) which is saturated withoxygen. The pulse oximeter consists of a probe attached to the patient'sfinger or ear lobe which is linked to a computerised unit.

SpO2 is typically reported as a percentage, values above 95% beingnormal and those below 95% indicating some degree of hypoxia (lack ofoxygen in the blood). Should a patient undergo an apnoea or hypopnoea,it is usual for the SpO2 signal to fall concomitantly with theventilation, albeit after some delay. During Cheyne-Stokes breathing theSpO2 signal will undergo the classic waxing and waning pattern alsocharacteristic of the ventilation.

Hence, it is conceivable that the algorithm described herein might use aflow signal estimate (ventilation) or an SpO2 signal or both signals toclassify breathing patterns as being typical of Cheyne-Stokes, OSA,mixed apnoeas etc.

FIG. 12 shows a Cheyne-Stokes patient's nasal flow signal over about 15minutes. FIG. 13 shows the same patient's SpO2 signal (saturation) andventilation signal (low-pass filtered absolute value of flow). Thesignals have been normalized and shifted to display them in the samegraph. The same pattern recognition techniques may be applied to bothsignals. For example: segment the signal into hypopnoeas/hyperpnoeas;analyze the shape of the hypopnoeas; determine the cycle lengths andspace ratios; perform a spectrogram (average of absolute value of anumber of FFTs); determine peaks jn the spectrogram at the CS frequency;determine a morphologic feature in both signals such as the jumpfeature; perform a continuous wavelet transform on both signals and useridge finding techniques to follow any CS frequency component over time.

EXAMPLE 1

A set of data for testing the ability of flow data to be classified intoOSA and CS instances consisted of 90 patient studies of approximately 8hours each. For purposes of the test, both nasal pressure, flow and twoeffort signals (abdomen & thorax) were recorded, making a confirmingdiagnosis of the underlying disease possible. The set was divided into 3groups of 30 patients: OSAi(obstructive apnoea), CS and Mixed. The datawere further classified (initially) on a 30-minute time-bin basis. Thetime periods were classified into the following categories: No apnoeasor hypopnoeas (<5) within the time period; Primarily CS breathing(>90%); Primarily OSA (>90%); Primarily (>90%) apnoeas and hypopnoeas ofmixed type (i.e. having a central component followed by a number ofobstructed breaths); A combination of different events, typically briefperiods of CS or mixed apnoeas interspersed between OSA; Patient ismoving and the signal is of too low a quality to be useable.

Typically if CS disease is present, CS breathing will occur in largeblocks of at least 20-30 minutes. The data set contained very fewperiods of “pure” mixed apnoeas. Rather, the mixed group of 30 patientscontained periods of OSA, CS breathing or a mixed picture.

Feature Analysis

All features were analyzed by calculating distributions for thedifferent groups (OSA, Mixed, and CS). The distribution was normalizedby application of an appropriate function, for example FIG. 5 shows thedistribution of the normalized max jump in hyperpnoea signal betweenbeginning of hyperpnoea and time of peak flow after application oflog-to-base-e. The leftmost curves represent a “normal” or Gaussiandistribution. It can be seen that the application of the log functionhas normalized the distributions and, further, that this feature showsgood separation between the CS (left) and OSA (rightmost) groups.

Cluster Analysis

Both k-means and fuzzy k-means clustering techniques were utilized tovisualize feature separation power. The features were first averaged ona per-patient basis and then custer analysis was used to demonstrate anatural clustering into correct groups. FIG. 6 shows such an analysis.The Euclidean 2-norm distance from each of two cluster centers isplotted one against the other. The CS and OSA patients naturally fallinto two groups except for one CS patient. The Mixed patients fall intoone group or the other depending on the length of time spent during thenight in different breathing patterns. The separating diagonal in thefigure represents a naive classifier suitable for per-patient grouping.What such a classifier cannot do is find a short period of CS breathingfrom amongst an otherwise OSA-dominated night.

Feature Temporal Averaging

The training of a classifier using patterns assigned to individualevents is problematical. Temporal averaging was used to reduce theamount of calculation, while also (potentially) increasing statisticalpower. A 30-minute time-bin was chosen as a best first-guess. Aftertemporal averaging, a new set of per-time-bin patterns is created. Theraw features used (visible separation of groups) were: hypopnoea length;hyperpnoea length; 1^(st) Fourier shape feature; 2^(nd) Fourier shapefeature; and normalized max jump. The time-averaged 30-minute binfeatures tested were (std=standard deviation, meansq=mean of square ofvalues, sqrt=square root, shift=allows calculation of temporaldifference).

Classifier Training and Testing

Once the data had been processed and the “expert” diagnosis made, agroup of 1440 30-minute bins was available for classifier training (90patients×16 bins).

Classifier Selection

Numerous statistical methods exist for the training of a classifier fromn-dimensional data, e.g.: nearest neighbor, neural nets, clusteranalysis. However, because the data “appeared” linearly separable,Bayesian decision theory was used. This theory (which relies onunderlying normal probability density functions) uses the minimizationof the Bayes error to calculate a discriminant surface. Such a surfaceseparates the data into one of n categories (in this case 2). Bothlinear and quadratic discriminant functions were utilized. The formerseparates the data with a hyperplane in m dimensions (where m is thenumber of features) while the latter separates the data with ahyperquadric. A hyperplane discriminant is always preferred (assumingaccuracy of the same order), as it will tend to be well behaved in areasof minimal data coverage.

Over-optimistic Train and Test

The classifier was trained using the training set and then theclassifier was tested using the same data. This results inover-optimistic values of sensitivity and specificity, as one wouldintuitively expect. However, again this is an insightful process and onecan use a minimal features set (≦3 features) in order to visualize theresult. FIG. 7 shows an LD classifier (plane shows equi-probabilitysurface). FIG. 8 shows a QD classifier (quadric shows equi-probabilitysurface).

Results

During each test the accuracy, sensitivity and specificity were noted aswas the current features set (or group of feature sets) with the bestaccuracy. Estimates of accuracy, sensitivity and specificity resulted ofthe order of 91%, 91% and 96% respectively.

EXAMPLE 2

Flow Filtering

The flow is filtered first to remove unwanted and uninterestinghigh-frequency content. The filter used is a digital FIR (finite impulseresponse) filter designed using the Fourier method using a rectangularwindow. The filter has a pass-band from 0 to 0.9 Hz, a transition bandfrom 0.9 to 1.1 Hz and a stop band above 1.1 Hz. The number of terms inthe filter varies with sampling frequency. The flow signal is filteredby convolving the time series point-wise with a filter vector.

Ventilation Calculation

A long-term ventilation signal y_(long) is calculated using a simple(first order) low-pass filter applied to the flow signal. A timeconstant of 200 seconds is used (longer than the longest possible cycleof Cheyne-Stokes breathing). In order to measure ventilation (and notmean flow), the filter is applied to the square of the flow signal andthe square root is taken of the filter output. Next, a ten-secondventilation y₁₀ is calculated (a more “recent” measure). This measure iscreated by convolving the square of the flow signal with a 10-secondsquare wave with an area of one, i.e. a 10-second-long moving average,and then taking the square-root of the result. This filter will have afive second delay constant over the frequency range of interest. Forthis reason the signal is shifted left by five seconds so that it “linesup” with the original signal for timing purposes. FIG. 10 shows aCheyne-Stokes flow waveform (large amplitude rapid varying curve), thelong-term ventilation (low amplitude slowly varying curve) and theleft-shifted 10-second ventilation (moderately varying curve) for atypical patient.

Event Detection from Ventilation Signals

The 10-second ventilation signal is used to create low and highthresholds for detection of events (hypopnoea-hyperpnoea sequences). Thethresholds are:thresh_(low)=0.5×y _(long);thresh_(high)=0.75×y _(long);The timing of events is calculated using the following algorithm:

-   1. Find all points where y₁₀<thresh_(low).-   2. Find all contiguous sections of the above points. These are    provisional hypopnoeas.-   3. Find all points where y₁₀>thresh_(high).-   4. Iterate over all of the hypopnoeas identified in step 2. If no    points identified in step 3 (hyperpnoeas) fall between hypopnoea n    and hypopnoea n+1, then the hypopnoeas n & n+1 are joined together    (because no hyperpnoea has been identified between them) to form one    hypopnoea. Repeat for all iterations.-   5. The hypopnoeas are now confirmed. All inter-hypopnoea regions are    considered hyperpnoeas. Each hypopnoea-hyperpnoea “event”    constitutes one possible Cheyne-Stokes cycle. E.g. in FIG. 10 there    are five cycles shown.    Calculate Event Timings

Event timings are calculated for each event as follows:τ_(hypopnoea) =t(end_of_hypopnoea)−t(beginning_of_hypopnoea);τ_(cycle) =t(beginning_of_next_hypopnoea)−t(beginning_of_hypopnoea);τ_(hyperpnoea) =t _(cycle) −t _(hypopnoea);

Obviously the above events will include some unwanted “garbage”. Forexample, a one-hour-long period of normal breathing bracketed on eachside by Cheyne-Stokes breathing will look like a one-hour-longhyperpnoea! (y10 always greater than threshold). Hence, the followingsensible limits are applied to the events:

-   τ_(hypopnoea): minimum=10 seconds, maximum=100 seconds,-   τ_(cycle): minimum=15 seconds, maximum=250 seconds,-   τ_(hyperpnoea): minimum=5 seconds.

All events outside these limits are rejected and not processed. We nowhave event timings and the ability to extract parts of the flow waveformfor further analysis. For example, we can iterate over all the eventsand select out those parts of the flow waveform that correspond tohyperpnoeas. FIG. 11 is an example where we have selected out anhyperpnoea from the above sequence and plotted it separately. In allfurther processing it is the 1 Hz filtered flow signal that is used forfeature extraction.

While the invention has been described in connection with what arepresently considered to be the most practical and preferred embodiments,it is to be understood that the invention is not to be limited to thedisclosed embodiments, but on the contrary, is intended to cover variousmodifications and equivalent arrangements included within the spirit andscope of the invention.

1. A method for diagnosing the presence of sleep disorders comprisingthe steps of: recording a signal representative of respiration from apatient using a logging device which includes a data-acquisition systemand a memory, processing the respiratory signal either on-board by therecording device or offline using a computer, dividing the signal into nepochs of equal length, recording events consisting of an hypopnoeafollowed by an hyperpnoea, detecting for each event its beginning andendpoints, calculating event lengths, and processing each hyperpnoea toderive shape features; wherein the shape features are calculated usingsingular value decomposition of the hyperpnoea ventilation signal by thesteps of extracting the hyperpnoea from the respiratory signal, forminga ventilation signal from the respiratory signal, scaling theventilation signal to give a vector of values, calculating shape factorsfrom the product of the pseudoinverse of the matrix and the vector ofvalues; and the matrix of orthogonal functions is:$M_{Basis} = \begin{pmatrix}{\frac{1}{\sqrt{\pi}}{\sin\left( \frac{t}{2} \right)}} \\{\frac{1}{\sqrt{\pi}}{\cos\left( \frac{t}{2} \right)}} \\\frac{\left( {{3\pi\;{\sin(t)}} - {8\;{\cos\left( \frac{t}{2} \right)}}} \right)}{\sqrt{\pi\;\left( {{9\pi^{2}} - 64} \right)}} \\\frac{\left( {{3\pi\;{\cos(t)}} - {4\;{\sin\left( \frac{t}{2} \right)}}} \right)}{\sqrt{\pi\;\left( {{9\pi^{2}} - 16} \right)}}\end{pmatrix}$ where t is the time base over the hyperpnoea from 0 to 2pi.
 2. A method for diagnosing the presence of sleep disorderscomprising the steps of: recording a signal representative ofrespiration from a patient using a logging device which includes adata-acquisition system and a memory, processing the respiratory signaleither on-board by the recording device or off line using a computer,dividing the signal into n epochs of equal length, recording eventsconsisting of an hypopnoea followed by an a hyperpnoea, detecting foreach event its beginning and end points, calculating event lengths,processing each hyperpnoea to derive shape features, and deriving a jumpfeature for each hyperpnoea by the steps of: extracting the hyperpnoeafrom the respiratory signal, forming a ventilation signal from anabsolute value of the respiratory signal, and approximating an envelopeof the ventilation signal with a droopy peak-detector.
 3. The method fordiagnosing the presence of sleep disorders of claim 2, wherein thefollowing droopy peak-detector is used to approximate the envelope ofthe ventilation signal: e[1] = v[1] for  i = 2…  m if  v[i] ≥ e[i − 1]e[i] = v[i] else${e\lbrack i\rbrack} = {{e\left\lbrack {i - 1} \right\rbrack} + {\frac{1}{2.5f_{s}}\left( {{v\lbrack i\rbrack} - {e\left\lbrack {i - 1} \right\rbrack}} \right)}}$end where e[i] is the approximate envelope, f_(s) is a samplingfrequency and v[i] is the ventilation signal, wherein an envelope isinterpolated over a time base to give e[i] a maximum positive differencein the interpolated signal in an interval between the beginning of theenvelope and the point at which the envelope attains its maximum valueis found, and the maximum difference by a mean value of the ventilationsignal to give a jump feature is scaled.
 4. A method for diagnosing thepresence of sleep disorders comprising the steps of: recording a signalrepresentative of respiration from a patient using a logging devicewhich includes a data-acquisition system and a memory, processing therespiratory signal either on-board by the recording device or off-lineusing a computer, dividing the signal into n epochs of equal length,recording events consisting of an hypopnoea followed by an hyperpnoea,detecting for each event its beginning and end points, calculating eventlengths, processing each hyperpnoea to derive shape features, andcalculating an additional feature from the entire epoch signal from aspectrogram of the epoch signal by: calculating a mean of therespiratory signal and subtracting said mean from the respiratorysignal, chopping the subtracted signal into n slices, windowing eachslice using a Hanning window, applying a Fourier transform to eachwindowed slice yielding a complex vector result for each slice,determining an absolute value of each complex result so as to yield areal valued vector per slice, averaging said real valued vectors toyield one vector, taking a natural log of a subsequent vector,extracting values in a frequency range of 0 Hz to 0.075 Hz to form asub-vector, de-trending the sub-vector and determining its mean, andcalculating the feature as a maximum minus the mean of said de-trendedsub-vector.